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We  have  accomplished  the  following  milestones: 

(1)  Developed  an  integrated  CAD-3DAM-FEM  approach.  The  integration  of  all  these  individual  components  of  a  CAE  approach 
to  the  analysis  of  material  and  structures  is  enabled  by  the  use  of  the  Python  programming  language. 

(2)  Developed  a  method  for  a  rational  approach  to  parameter  calibration  in  TIM  model  systems.  We  have  defined  this  approach 
to  determine  the  contact  stiffness  between  unit  elements  as  well  as  the  friction  conditions.  Excellent  agreement  between  model 
and  simulation  is  obtained. 

(3)  Provided  understanding  of  the  high  rate  loading  (constant  rate)  of  TIM  assemblies.  We  have  found  that  the  principle 
underlying  mechanism  of  the  formation  of  force  chains  is  again  present,  albeit  somewhat  mitigated  if  slip  is  considered.  High 
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such  interaction  exist. 

(5)  Investigate  the  influence  of  local  material  failure  on  TIM  failure.  We  embed  cohesive  zone  models  for  the  description  of 
material  failure  into  the  analysis.  We  find  that  residual  velocities  are  potentially  favorably  affected  by  fragmentation.  We  also 
determine  a  size  effect  on  the  predicted  residual  velocities.  TIM  systems  response  is  found  to  be  dependent  on  the  material 
internal  length  scale  as  introduced  by  the  cohesive  zone  model. 
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1.  Introduction 


In  the  design  of  multifunctional  structural  components  it  is  essential  to  go  beyond  the 
selection  of  materials  with  the  appropriate  properties.  Developing  logical  and  hybrid 
structures  can  be  effective  in  solving  some  of  critical  and  challenging  engineering  problems 
that  may  require  multifunctionality  without  developing  new  and  exotic  material  composition. 
Hybrid  material  systems  are  understood  to  provide  efficient  alternatives  to  conventional 
materials  and  fill  holes  in  material  property  spaces  [1-3].  One  of  the  methods  relevant  for  the 
creation  of  hybrid  material  systems  is  that  of  segmentation.  Thereby,  a  hybrid  material 
system  is  created  from  the  bottom-up  by  the  ordered  assembly  of  unit  elements.  Hybrids 
created  by  this  approach  allow  for  the  exploration  of  important  material  mechanics  concepts 
such  as  topological  toughening  and  of  size  effects.  Past  studies  on  hybrid  materials  created  by 
assembly  have  predominately  considered  quasi-static  loading  conditions.  The  aim  of  the 
present  study  is  to  explore  the  response  of  such  material  systems  to  high  rate  loading  and  to 
impact.  In  particular,  the  study  is  concerned  with  one  class  of  segmented  hybrids, 
Topologically  Interlocked  Material  assemblies,  TIMs.  In  a  TIM  system,  each  individual  unit 
element  is  topologically  confined  by  its  neighbors  due  to  the  oblique  angles  of  the  basic  unit 
elements;  therefore  no  binders  are  necessary  to  maintain  the  integrity  of  the  assembly  [4-8], 
We  hypothesize  that  the  topological  interlocking  of  the  unit  (or  fragment)  elements  would 
provide  for  an  enhancement  in  the  impact  resistance  compared  to  monolithic  solids. 

The  concept  of  topological  interlocking  assemblies  has  attracted  some  interest  in  the  recent 
past.  Proposed  by  Glickman  [8]  with  the  name  of  G-Block  for  the  reduction  of  sub-structure 
and  diminished  maintenance  cost  for  paving  system,  expanded,  generalized  and  termed  with 
the  current  name  by  Dyskin  et  al.  [4,  9,  10],  TIMs  can  be  considered  as  novel  mechanical 
meta-materials  as  these  derive  their  unique  properties  not  from  composition  but  structure. 
Forces  and  energy  can  only  be  transmitted  through  the  contact  surfaces;  accordingly,  there 
are  no  tensile  forces  developed  when  TIMs  carry  loads.  These  unique  characteristics  of  TIMs 
provide  unique  properties  that  can  be  utilize  to  create  adaptive  and  configurable  structures  to 
harsh  conditions  such  as  random  and  harmonic  vibrations,  thermal  loads,  repetitive  shocks 
and  acoustic  attenuation.  For  example,  Dyskin  et  al.  [4,  5,  7,  9,  10]  introduced  the  general 
concept  of  TIM,  and  demonstrated  the  requirements  for  basic  elements  to  assemble  a  TIM  tile. 
Estrin  et  al  .  [11]  assembled  TIMs  with  identical  cubes  and  tested  their  responses  with  point 
load  tests.  Schaare  et  al.  [12]  investigated  the  response  of  cube-based  TIMs  to  point  loading 
both  experimentally  and  numerically.  Brugger  et  al.  [13]  conducted  indention  experiments  on 
TIMs  assembled  by  osteomorphic  ice  blocks  and  plaster  made  cubes.  Their  results  showed 
that  TIM  composed  of  cube  elements  can  experience  negtive  stiffness  when  unloading. 
Dyskin  et  al.  [14]  revealed  that  TIMs  are  damage  tolerant  in  that  cracks  inside  of  one  unit 
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element  will  not  propagate  across  the  contact  interface  into  adjacent  units.  Khandelwal  et  al. 
[15]  experimentally  examined  the  quasi-static  responses  of  TIMs  asembled  with  regular 
tetrahedra,  and  proposed  an  analytic  model  to  predicted  the  quasi-static  beaviors  of  TIM  with 
thrust  line  theory.  Mather  et  al.  [16]  studied  the  remanufacturability  of  TIMs,  illustrated  that 
only  small  portion  of  units  could  not  be  reused  after  failure  of  the  TIMs.  And  for  applications 
of  TIMs,  Estrin  et  al.  [17]  proposed  the  potential  application  of  TIM  as  protective  tiles  for  the 
space  shuttle.  Carlesso  et  al.  [18]  illustrated  TIM  composed  by  porous  osteomorphic  blocks 
enhances  sound  absorption  because  of  the  gaps  among  basic  blocks. 

The  present  investigation  is  aimed  at  understanding  the  response  mechanism  of  TIMs  due  to 
impact  loading,  and  the  investigation  of  potential  benefits  of  TIMS  in  energy  absorption,  as 
well  as  mechanical  and  structural  applications.  TIMs  consists  of  dense  packing  of  tetrahedral 
unit  elements  made  of  Acrylonitrile  Butadiene  Styrene  (ABS)  material,  which  were 
fabricated  using  fused  deposition  modeling  (FDM)  additive  manufacturing  (AM)  motivate 
the  material  system  under  investigation  [15].  Limited  exploratory  research  has  been  done  in 
understanding  the  mechanical  characteristics  of  these  novel  structures,  more  specifically 
TIMs  subject  to  low  velocity  impact.  Plates  subjected  to  impact  are  of  great  research  interest 
for  protection  purpose,  and  extensive  work  has  been  reported  on  different  types  of  systems 
subjected  to  ballistic  impact  [19-27].  Most  recently,  Wadley  et  al.  [28]  studied  sandwich 
structures  with  different  shapes  of  core  material  under  ballistic  impact.  Their  results 
demonstrate  that  sandwich  structure  performance  can  be  improved  with  topological 
optimization  of  core  materials.  In  the  present  study,  the  impact  velocities  considered  are 
much  lower  than  the  dilatational  wave  speed  of  the  material  in  the  unit  elements  (V!mp  < 
Cd/10).  The  finite  element  method  was  employed  to  conduct  numerical  experiments.  First,  the 
model  is  calibrated  on  experimental  data  of  low  and  constant  velocity  loading.  Then,  the 
response  under  constant  applied  velocity  is  considered,  and  the  contributions  of  inertia, 
elastic  resistance,  and  contact  and  friction  to  the  mechanical  resistance  are  explored. 
Subsequently,  impact  loading  is  considered  and  the  formulations  of  Recht  and  Ipson  [29]  as 
well  as  Lambert  and  Jonas  [30]  are  employed  to  characterize  the  velocity  response.  These 
models  are  reformulated  to  provide  specific  insight  into  the  impact  response  of  the  TIM 
assemblies.  A  comparison  to  typical  data  on  monolithic  structures  under  impact  is  made. 

2.  Python  Scripting  Based  Model  Creation  and  Manufacturing 

The  ABAQUS  Scripting  Interface  enables  the  user  to  effectively  customize  the  creation  of 
finite  element  models.  With  the  API,  one  can  automate  repetitive  tasks,  extend  functionality 
and  enhance  the  interface  through  the  use  of  Python.  What  makes  this  scripting  language  so 
powerful  is  the  fact  that  for  every  feature  in  ABAQUS /CAE,  there  is  a  corresponding  Python 
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script  command  that  can  be  modified  to  provide  greater  functionality  and  more  effective 
automation.  In  the  present  project,  Python  scripts  were  used  to  create  ABAQUS/CAE  models, 
define  material  properties,  generate  meshing,  apply  boundary  conditions  and  loads,  and  to 
submit  job  for  analysis.  Furthermore,  the  Python  scripts  were  also  employed  to  generate  CAD 
models  for  additive  manufacturing.  Granular  crystals  are  considered  as  assemblies  of  many 
identical  unit  elements.  Therefore,  many  repetitions  in  the  model  creation  would  be  needed  if 
the  process  was  conducted  without  a  script.  The  granular  crystal  considered  here  was  an 
assembly  of  interlocking  tetrahedra  arranged  in  a  dense  packing  in  the  plane  [31].  The 
granular  crystal  -  also  called  a  Topologically  Interlocking  Material  (TIM)  [10]  -  is  confined 
at  its  boundaries  by  rigid  abutments,  and  a  projectile  type  loading  is  applied  at  the  center  of 
the  assembly  via  a  rigid  body  (see  Figure  1). 


Figure  1.  CAD  model  of  the  TIM  plate  and  the  projectile 
The  general  procedure  for  generating  simulation  with  Python  script  is  as  follows: 

(1) .  Create  a  model  database  named  as  tetrahedron 

(2) .  Create  a  part  named  as  tet,  which  is  an  individual  tetrahedron  element  of  the  TIM  plate. 

Define  the  material  properties  for  the  tetrahedron  and  assign  them  to  the  part  of  tet.  Set 
seeds,  mesh  the  tetrahedron  part,  assemble  the  part,  pattern  the  instance  and  position 
them  by  translation  and  rotation. 

(3) .  Create  a  part  named  as  boundary,  which  is  one  side  of  the  support  frame.  Define  the 

material  properties  and  assign  them  to  the  part  of  boundaiy.  Set  seeds,  mesh  the  frame 
part,  assemble  the  part,  pattern  the  instance  and  position  them  by  translation  and  rotation. 
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(4) .  Create  a  part  named  as  indent er,  which  is  used  to  impact  the  TIM  plate.  Define  the 

material  properties  and  assign  them  to  the  part  of  indenter.  Set  seeds,  mesh  the 
tetrahedron  part,  assemble  the  part,  and  position  it  by  translation  and  rotation. 

(5) .  Define  the  contact  interactions  and  rigid  bodies.  In  the  simulation,  the  indenter/projectile 

and  the  support  frame  are  considered  as  rigid  bodies. 

(6) .  Apply  boundary  conditions  and  loads.  For  the  support  frame,  all  degrees  of  freedom  are 

constrained.  A  predefined  velocity  is  applied  on  the  projectile,  which  will  impact  the 
TIM  plate  at  the  plate  center. 

(7) .  Submit  the  job  for  analysis. 

Parameterization  is  a  distinct  advantage  of  Python  scripting  for  ABAQUS  as  it  allows  for 
high  efficiency  in  modeling,  particularly  for  dimensional  analysis.  For  example,  Figure  2 
shows  a  TIM  tile  with  7x7  tetrahedra.  To  have  a  TIM  tile  with  10x10  tetrahedra  (Figure  3), 
we  just  need  to  set  the  number  of  tetrahedral  in  one  row/column  to  be  10. 


Figure  2:  TIM  plate  with  7x7tetrahedra 

Two  basic  dimensions  for  the  TIM  plate  are  the  edge  length  of  individual  tetrahedron  and  the 
numbers  of  tetrahedra  in  one  row/column.  Therefore  these  two  variables  are  defined  as  the 
basic  dimensional  parameters  for  the  analysis  in  the  Python  scripts  and  all  other  dimensions 
can  be  obtained  based  on  them.  All  the  material  properties  are  also  parameterized,  which 
means  that  the  material  properties  can  be  changed  very  easily. 
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Figure  3:  TIM  plate  with  10x10  tetrahedra 

Cracks  may  initiate  and  propagate  in  materials  when  they  are  subjected  severe  loading 
conditions.  Many  approaches  have  been  developed  to  predict  the  process  of  crack  initiation 
and  propagation  in  the  frame  of  fracture  mechanics  [32,  33].  Abaqus  has  implemented 
cohesive  zone  method  (CZM),  virtual  crack  closure  technique  (VCCT),  and  extended  finite 
element  method  (XFEM)  to  simulate  the  crack  propagation.  With  the  current  release  of 
Abaqus  package  (Abaqus  6. 12),  all  the  three  methods  can  work  under  Abaqus/standard,  while 
only  the  CZM  can  be  used  under  Abaqus/explicit.  Therefore  CZM  is  employed  for  the 
current  project  to  predict  the  damage  of  topologically  interlocked  material  (TIM)  subjected  to 
impact.  Moreover,  the  microstructure  nature  of  the  tetrahedra  to  assemble  the  TIM  tile,  which 
are  manufactured  by  3D  printing,  is  very  desirable  to  the  CZM. 

On  the  other  hand,  the  creation  of  CZM  may  be  very  tedious  if  possible,  particularly  when  the 
models  are  complex,  since  the  modification  of  nodal  coordinates  is  necessary  to  have  zero 
thickness  cohesive  layers.  With  its  powerful  loop,  condition  and  selection  statements,  Python 
scripting  can  dramatically  improve  the  efficiency  of  modeling  CZ  with  Abaqus.  Furthermore, 
with  the  help  of  parameterization,  it  is  very  convenient  to  modify  the  number  of  CZ  layers  in 
a  tetrahedron,  or  even  the  orientation  of  the  CZ  layer. 

The  additional  steps  for  CZ  enhanced  model  are  to  define  the  CZ  properties  and  create 
cohesive  elements,  comparing  with  the  models  without  cohesive  connections.  In  the 
simulations,  the  bilinear  traction-separation  law  is  used  to  describe  the  CZ  constitutive 
behavior.  In  choosing  the  CZ  properties,  such  as  failure  displacement,  stiffness  and  softening, 
both  prediction  accuracy  and  computational  resources  need  to  take  into  consideration  [34]. 

Figure  4  illustrates  the  procedure  of  creating  the  cohesive  elements  in  a  tetrahedron  part.  First, 
layered  tetrahedron  part  is  created  (Figure  4a)  and  meshed  (Figure  4b).  To  have  zero 
thickness  cohesive  zones,  the  cohesive  element  nodes  are  modified  (Figure  4c). 
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Figure  4:  Cohesive  zone  model  of  tetrahedron 

(a)  Geometry  model;  (b)  mesh  without  modification;  (c)  mesh  with  modification 

Python  scripts  can  also  be  applied  to  the  creation  of  CAD  model  and  STL  files  for  the  3D 
print  manufacturing  of  granular  crystals.  In  this  application  the  Python  script  creates  ACIS 
SAT  3D  model  files  (.sat),  which  then  further  are  saved  in  the  STL  file  format  for  processing 
with  3D  printers  such  as  the  Objet  Connex  350  (Stratasys).  Figure  5  shows  the  CAD  model 
for  the  TIM  plate  with  7x7  tetrahedra.  Figure  6  shows  the  importation  to  CAD  software  of 
the  CAD  model  created  by  the  Python  scripts.  Figure  7  a  screen  shot  of  the  printing  machine 
interface  software,  and  Figure  8  the  final  3D  printed  prototype. 


Figure  5:  CAD  model  for  the  TIM  plate  with  7x7  tetrahedra 
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Figure  6:  Importing  the  CAD  file  to  CAD  software 
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Figure  7:  Printing  the  CAD  model  with  printing  machine  interface  software 


Figure  8:  3D  printed  prototype  of  a  TIM  plate  with  7x7tetrahedra 

3.  FEM  Model  Definition  and  Calibration 

3.1.  Model  Definition 

Figure  9  shows  the  model  system  under  consideration  with  (a)  a  regular  tetrahedron  as  the 
basic  unit  of  the  TIM  assembly,  (b)  a  boundary  element  (termed  abutment),  (c)  the 
indenter/projectile  and  the  center  five  tetrahedra  onto  which  the  indenter  transfers  the  load,  (d) 
the  TIM  tile  composed  with  identical  regular  tetrahedra,  (e)  the  overall  TIM  assembly 
together  with  the  cylinder  shaped  indenter/projectile.  In  the  model  assembly,  total  of  N><N=  49 
identical  tetrahedra  are  assembled  in  a  square  lattice  pattern  into  a  monolayer  TIM  assembly 
[4,  8,  15], and  confined  with  abutments  along  the  edges  of  the  assembly.  The  edge  length  of 
an  individual  regular  tetrahedron  is  so.  The  other  dimensions  of  the  TIM  assembly  can  be 
derived  as  the  thickness  of  the  monolayer/  =  ,s'o/V2,  the  dimension  of  the  mid-plane  of  the  TIM 
assembly ,L=sdx(N/2).  Here,  the  TIM  assembly  is  characterized  by  N  =  7  andvo  =  25  mm, 
such  that  t  =  17.7  mm  and L  =  87.5  mm. 


Figure  9.  Model  of  Topologically  Interlocked  Material  (TIM)  assembly:  (a)  Unit  element 
regular  tetrahedron,  (b)  Abutment  (c)  Indenter/projectile  and  center  tetrahedra,  (d)  TIM  tile, 

(e)  Overall  assembly. 

Two  loading  conditions  were  considered:  (1)  constant  velocity  loading  by  an  indenter,  and  (2) 
impact  loading  by  a  projectile.  In  both  cases  the  load  was  applied  to  the  center  of  the  TIM 
assembly  such  that  the  indenter/projectile  interacts  directly  with  tetrahedra  #1,  #2,  and  #3  of 
the  assembly.  The  diameter  D  of  the  indenter/projectile  is  D  =  L  /  5  . 

Tetrahedra  were  considered  to  be  isotropic  and  elastic.  Motivated  by  prior  experiments  on 
TIM  assemblies  in  [15],  the  material  properties  of  ABS  were  employed  (Young's  modulus 
£'=1.827  GPa,  Poisson's  ratio  v  =  0.35,  density  p0  =  950  kg/m3).  One  tetrahedron  thus 
possesses  a  mass  of  mo=l. 75x10°  kg  and  the  overall  assembly  weighs  8.58xl0"2  kg.  The 
abutments  as  well  as  the  indenter/projectile  were  considered  as  rigid  bodies.  The  density  of 

the  projectile  was  p  =1470  kg/m3.  The  contact  properties,  i.e.,  the  contact 

pressure-overclosure  relation  and  the  friction  condition  were  not  directly  available.  A 
Coulomb  friction  model  was  employed  which  introduces  the  coefficient  of  friction  p.  For  the 
normal  contact  behavior,  it  was  assumed  that  contact  pressure  pc  is  linearly  dependent  on  the 
overclosure  An  as/;c  =  Kc Am  where  Kc  =  (K*E*)/ (ns  o)is  the  contact  stiffness  with  K*  a 
normalized  contact  stiffness,  and  E*  =  £/[2(  l  -v2)].  This  type  of  interactions  was  adopted  to 
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describe  the  contact  behaviors  for  the  inter-tetrahedra,  tetrahedron-abutment  and 
tetrahedra-ind enter  contacts.  Experimental  data  on  the  low  velocity  impact  response  of  TIM 
assemblies  presented  in  [15]  was  used  to  calibrate  the  model  parameters  jl  and  K  . 

In  the  constant  velocity  cases,  a  predetermined  displacement  is  applied  on  the  indenter  and 
the  time  interval  of  loading  is  varied  in  order  to  achieve  a  constant  velocity  Vo.  In  the  impact 
cases,  a  projectile  of  a  defined  mass  Mwith  initial  velocities  Vunp  strikes  the  TIM  assembly. 
The  mass  of  the  projectile  is  chosen  as  10%  of  the  total  mass  of  the  7X7  tetrahedra.  The 
abutments  are  fixed  in  space  in  all  cases. 

Finite  element  models  of  the  TIM  assembly  were  constructed  in  the  general  purpose  finite 
element  code  ABAQUS  (Version  6.12)  [35].  The  finite  element  mesh  employs,  3D  linear  solid 
element  (C3D8R  in  the  ABAQUS  code)  with  each  tetrahedron  is  comprised  of  500  elements 
such  that  the  overall  model  consisted  of  24,500  elements.  The  model  was  solved  using  the 
ABAQUS/Explicit  solver,  and  the  general  contact  algorithm  implemented  in  the  ABAQUS 
software  was  employed.  In  all  the  simulations,  one  analysis  step  was  defined.  The  automatic 
time  increment  algorithm  was  employed  and  the  global  stable  increment  estimator  was 
selected.  The  step  time  period  for  constant  velocity  loading  varies  dependent  on  loading 
velocities  while  for  impact  loading  an  analysis  interval  of  5  ms  is  considered  in  all  the  cases. 

3.2.  Parameter  Calibration 

In  the  experiments  presented  in  [15],  a  drop  mass  of  Mdrop=  6.21  kg  and  drop  height  Hdr0p= 
36.58  mm  were  employed.  Since  the  drop  mass  significantly  exceeds  that  of  the  TIM 
assembly,  a  near  constant  velocity  loading  condition  (Vo  =  0.925  m/s)  was  achieved  in  the 
experiment.  The  data  of  the  forces  -  deflection  record  were  presented  in  [15],  and  this  data 
was  used  to  calibrate  the  parameters  jl  and  K*. 

The  computational  study  provided  the  force-deflection  curves  of  TIM  for  contact  stiffness  ( K 
=  0.5)  and  different  coefficients  of  friction.  Figure  10  shows  predicted  force-deflection  curves 
from  computations  considering  a  range  of  values  of  jl  (0.1<  jt  <1.0).  In  these  calculations 
the  contact  stiffness  was  K*  =  0.5.  In  all  cases  a  force  -  deflection  response  typical  of  TIM 
assemblies  was  predicted  [15].  The  force  initially  increases  with  the  deflection;  a  peak  force 
value  is  reached,  followed  by  gradual  softening  until  the  resistance  of  the  TIM  assembly 
vanishes.  The  resistance  of  the  TIM  assembly  to  the  applied  load  was  found  to  increase  with 
an  increase  in  jl.  The  predicted  peak  force  increased  from  48.5  N  to  513.0  N  when  jl  was 
changed  from  jl  =  0. 1  to  1.0.  The  results  also  indicate  an  increase  of  the  deflection  to  final 
collapse  of  the  TIM  with  an  increase  in  jl.  For  the  lowest  value  of  ^considered  (^=0.1),  the 
load  carrying  capacity  of  the  TIM  vanishes  at  a  deflection  of  22.0  mm  but  this  limiting 
deflection  more  than  doubles  for  cases  considering  higher  values  of  jl.  The  next  set  of 
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computations  considered  a  range  of  contact  stiffness  values  and  constant  coefficient  of 
friction.  Figure  1 1  shows  reaction  the  predicted  force-deflection  records  for  computations 
considering  variations  in  contact  stiffness  (0.1<if*<1.0).  In  these  calculations  //  =  0.3  was 
considered.  It  is  found  that  the  magnitude  of  the  predicted  forces  increases  with  increased 
contact  stiffness.  The  peak  force  was  found  to  increase  from  approximately  52.0  N  to  5 13.0  N 
as  K*  was  increased  from  K*=  0.1  to  1.0.  However,  the  deflection  to  loss  of  load  carrying 
capacity  did  not  vary  significantly,  and  for  all  the  cases  was  about  35.0  mm. 

The  results  of  the  parametric  computations  indicate  that  the  coefficient  of  friction  affects  both 
the  magnitude  of  the  force  and  the  deflection  to  final  failure,  while  the  contact  stiffness  plays 
a  role  only  in  the  force  magnitude.  Therefore,  the  following  calibration  procedure  was 
followed.  First,  the  coefficient  of  friction  was  calibrated  with  respect  to  the  experimental 
results  of  the  deflection  to  final  failure.  Then,  the  contact  stiffness  is  selected  to  provide  a 
final  force  -  deflection  response  which  fits  that  of  the  experiments  presented  in  [15].  The 
calibrated  parameters  are  //  =0.2  and  K*  =  0.38.  Model  predictions  of  the  force  -  deflection 
response  obtained  using  the  best-fit  parameters  for  ju  and  K*  together  with  the  experimental 
data  are  shown  in  Figure  12,  and  a  good  quantitative  agreement  is  found.  The  residual 
deformations  of  the  TIM  assemblies  after  indention  are  shown  for  the  experiment  in  Figure 
13(a)  [15],  and  the  model  simulation  in  Figure  13(b).  A  good  qualitative  agreement  is 
obtained. 


Figure  10.  Predicted  force-deflection  curves  for  a  range  of  values  of  the  coefficient  of  friction 


li 


Figure  1 L  Predicted  force-deflection  curves  for  a  range  nonnalized  contact  stiffness  values 

K\ 


Figure  12.  Force  -  deflection  curves:  experimental  data  [15]  and  best  fit  prediction. 
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(b) 

Figure  13.  Deformed  configurations  of  TIM  assemblies  after  loading:  (a)  Experimental 
observation[15],  (b)  Model  prediction  with  color  contours  depicting  the  residual  stress  state 
as  characterized  by  the  von  Mises  equivalent  stress. 

4.  Dynamic  Responses  of  TIM 

4.1.  Constant  Velocity 

Simulations  of  TIM  subjected  to  constant  velocity  loading  were  conducted  with  applied 
velocities  ranging  from  0.1  m/s  to  10  m/s.  Figure  14  depicts  the  predicted  force-deflection 
responses.  For  the  lowest  velocity  case,  Fo=0.1  m/s,  the  force  increases  smoothly  with  the 
increase  of  deflection,  reaches  its  peak  at  a  deflection  of  1 1.0  mm  then  decreases  gradually  to 
zero.  As  the  applied  velocity  increases,  multiple  peaks  develop  in  the  force-deflection  curves. 
When  the  loading  speed  is  greater  than  1  m/s  but  less  than  5  m/s,  the  first  peak  increases  with 
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increasing  velocity  dramatically,  while  the  second  one  maintains.  When  the  loading  speed  is 
greater  than  5  m/s,  both  peaks  increase  with  increasing  velocity. 

Figure  15  shows  the  force-deflection  relations  for  Vo  =  0.1  m/s  with  confinement,  Vo  =  4  m/s 
with  and  without  confinement,  respectively.  A  simulation  for  unconfined  TIM  means  the 
abutments,  shown  in  Figure  9(d),  are  taken  away  in  the  computation.  In  this  case,  the  force 
increase  quickly  to  its  peak  then  drops  precipitously.  There  is  no  support  force  acting  on  the 
TIM  tile,  therefore,  the  force  shown  in  the  figure  is  due  to  the  ’inertia'  of  the  TIM  tile.  Since 
there  is  no  such  phenomenon  occurred  for  the  Vo  =  0.1  m/s  case,  the  Vo  =  0.1  m/s  case  can  be 
considered  as  quasistatic  loading  case.  Therefore,  as  shown  in  Figure  15,  the  force-deflection 
curve  for  the  Vo  =  4  m/s  with  confinement  can  be  considered  as  the  sum  of  the  other  two 
cases:  quasi-static  response  and  unconfined  model  with  same  impact  velocity.  Details  of  the 
determination  of  the  inertia  of  TIM  is  discussed  in  Appendix. 
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— VQ  =  0.1  m/s  (Quasistatic) 


Figure  15.  Predicted  force  -  deflection  curves  for  quasistatic  loading  with  confined  model, 

Vo=  4.0  m/s  confined  and  unconfined  models. 

Figure  16  depicts  the  velocities  of  center  tetrahedra  during  the  indentation  for  relatively  low 
and  high  indention  speeds,  Vo  =  4  and  10  m/s,  respectively.  The  velocity  of  center  tetrahedron 
is  defined  on  the  center  of  the  center  tetrahedron,  and  the  center  tetrahedron  is  the  tetrahedron 
that  the  indenter  directly  strikes  on,  T1  as  shown  in  Figure  9(c).  The  velocities  increase  with 
time  at  the  beginning  of  indention,  then  decrease  and  increase  alternatively,  and  finally  they 
keep  constant.  Overshoots  occur  during  indention.  For  F0=  1 0  m/s,  after  its  peak  value,  the 
velocity  of  center  tetrahedron  drops  far  below  the  indentation  speed  because  of  the  interaction 
with  adjacent  tetrahedra  through  contact.  The  indenter  accelerates  it  again,  which  results  in  a 
second  peak  force  in  the  force-deflection  curves.  For  Vo  =  4  m/s,  there  also  has  time  that  the 
center  tetrahedron  speed  is  lower  than  that  of  indentation  speed.  However,  the  speed 
difference  is  not  significant  so  the  re-acceleration  is  not  observed  in  the  force-deflection 
curve. 
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Figure  16.  Predicted  translational  velocities  (in  load  direction)  of  the  center  tetrahedron  for 

loading  speeds  of  Vo=4  m/s  and  10  m/s. 

4.2.  Impact 

All  series  of  computations  with  impact  velocities  ranging  from  1  m/s  to  80  m/s  was 
performed,  all  with  employing  the  calibrated  model  parameters.  The  ballistic  limit  for  the 
TIM  assembly  is  determined  as  Vb/=  34  m/s,  such  that  at  the  end  of  this  impact  event  the 
velocity  of  the  projectile  is  zero.  Below  Vbi  the  projectile  rebounds  and  above  F/,/the  target  is 
penetrated.  Figure  17  (a)  and  (b)  show  the  force  -  displacement  responses  of  the  projectile  for 
low  and  relatively  high  impact  velocities,  respectively.  The  forces  increase  at  the  beginning 
of  impact,  reach  their  peaks  then  decrease  with  displacement.  Similar  to  the  constant  velocity 
loading  cases,  the  reaction  force  increases  with  increasing  of  loading  speed,  and  second  peaks 
developed.  While  for  the  same  initial  loading  velocities  (Figure  14  and  Figure  17(a)),  the 
impact  loading  results  in  lower  reaction  force.  For  higher  impact  cases  (Figure  17(b)), 
rebound  and  penetration  can  be  observed  for  different  impact  velocities.  Figure  17  also 
indicates  that  the  TIM  absorbs  higher  amount  of  impact  energy  with  higher  impact  velocity. 
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Figure  17.  Force  -  displacement  response  of  projectile  with  different  impact  velocities:  (a) 
Low  impact  velocity,  (b)  Relatively  high  velocity. 


Figure  18.  Predicted  translational  velocities  (in  load  direction)  of  the  center  tetrahedron  and 
projectile  for  loading  speeds  of  Vimp=  4  m/s  and  10  m/s. 

Figure  18  predicts  the  velocities  of  center  tetrahedron  and  projectile  for  impact  velocity  Vimp 
=  4  and  10  m/s,  respectively.  Comparison  of  the  results  with  constant  velocity  loading 
(Figure  16)  shows  that  center  tetrahedron  obtains  smaller  velocity  under  impact  loading  when 
the  initial  loading  speed  for  both  cases  are  the  same;  and  after  the  center  tetrahedron 
decelerated  due  to  interaction  with  surrounding  tetrahedra,  the  velocity  difference  between 
projectile  and  center  tetrahedron  is  also  smaller  in  the  impact  loading  case,  which  interprets 
the  differences  among  force  peaks  shown  in  Figure  14  and  Figure  17(a). 

Figure  19  depicts  the  velocity  histories  of  the  projectile  and  the  center  five  tetrahedra,  Figure 
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9(c),  for  four  values  of  the  impact  velocity  Vimp.  In  all  cases  Vimp>  Vti.  At  the  ballistic  limit, 
Figure  19(a),  both  velocities  of  the  tetrahedra  and  the  projectile  reach  zero.  For  impact 
velocity  slightly  higher  than  ballistic  limit,  Vimp=  45  m/s,  the  velocities  of  the  center 
tetrahedron  and  the  projectile  velocity  are  only  slightly  larger  than  that  of  the  surrounding 
tetrahedra,  Figure  19(b).  This  indicates  that  the  projectile  strongly  interacts  not  only  with  the 
center  tetrahedron,  but  due  to  the  topological  interlocking,  also  with  the  surrounding.  For 
relatively  high  impact  velocity  cases  shown  in  Figure  19(c)  and  (d)  this  situation  changes.  In 
both  cases  shown,  i.e.,  Vimp  =  55  and  80  m/s,  the  center  tetrahedron  gains  much  higher 
velocity  than  its  neighbors.  Still,  Vimp  =  55  m/s  the  projectile  velocity  is  similar  to  that  of  the 
surrounding  unit  elements,  indicating  that  the  TIM  assembly  continues  to  provide  some  level 
of  resistance.  Only  at  much  higher  impact  velocity,  Vimp  =  80  m/s,  is  such  an  interaction  lost 
and  the  projectile  velocity  increases  beyond  that  of  the  surrounding  tetrahedra.  The  images 
of  the  defonned  configuration  of  the  TIM  assemblies  during  impact,  Figure  20,  confirm  these 
different  levels  of  interaction  between  projectile  and  the  TIM  target. 


Figure  19.  Predicted  projectile  velocity  during  impact  compared  to  the  velocities  of  the  five 
central  tetrahedra  units  in  the  TIM  assembly:  (a)  Vimp  =  34m/s  (b)  Vimp  =  45  m/s  (c)  Vimp  =  55 

m/s  (d)  Vimp  =  80  m/s. 
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Figure  20.  Predicted  deformed  configurations  of  the  TIM  assemblies  at  time  =  5  ms: 

(a)  Vimp  =  34  m/s  (b)  Vimp  =  45  m/s  (c)  Vimp  =  55  m/s  (d)  Vimp  =  80  m/s. 

The  central  variable  to  characterize  the  resistances  of  plates  and  plate-like  systems  to  impact, 
is  the  residual  velocity  of  the  projectile  after  passing  the  target.  When  the  impact  velocity  is 
greater  than  the  ballistic  limit  velocity,  the  projectile  penetrates  the  target  and  keeps  moving 
forward  with  non-zero  certain  kinetic  energy.  For  conventional  metal  plates  or  plate-like 
targets, 

Lambert  and  Jonas  [30]  proposed  a  fonnula  to  predict  the  residual  velocity: 


'0,  Vimp  <  vbl 

a(vp  -vp)  V  >  V 

u  V  imp  v  bl  )  v  imp  ^  v  bl 


(0) 


in  which  a  an  p  are  constants  to  be  fitted  with  the  residual  velocity  data,  Vimi,  is  impact 
velocity  and  Vbi  is  ballistic  velocity.  Figure  21  depicts  the  fit  to  the  residual  velocity  data 
using  the  Lambert-Jonas  models.  A  nonlinear  least  square  method  was  employed  to  obtain 
the  parameters,  which  was  conducted  with  MATLAB  built-in  function  lsqnonlin  [36] .  The 
fitted  parameters  are  a  =  0.343  and  p  =  1.487.  The  coefficient  of  determination  (R2)  is  0.9570, 
indicating  the  curve  well  fits  the  numerical  results. 
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30 


•  FEA  results 


Figure  2L  Predicted  residual  velocities  in  dependence  of  the  impact  velocity  and  data  fit  by 

the  Lambert- Jonas  model. 


The  fitted  Lambert-Jonas  model  can  provide  an  overall  good  description  of  the  numerical 
results.  Still,  but  for  higher  impact  velocities,  the  model  fits  a  less  accurate  than  for  velocities 
close  to  the  ballistic  limit.  A  detailed  examination  of  Figure  21  leads  to  the  conclusion  that 
the  residual  velocities  solution  can  be  separated  into  at  least  two  regimes  (below  and  above 
56  m/s).  Such  differences  in  the  response  can  be  traced  back  to  the  velocity  history  of 
projectile  and  tetrahedra,  Figure  19,  and  is  due  to  the  differences  in  interaction  with  the  unit 
elements  in  the  TIM  assembly.  In  summary,  a  two  limiting  conditions  can  be  defined:  (i)  low 
velocity  regime  with  strong  persisting  projectile  -  target  interaction,  and  (ii)  a  high  velocity 
regime  with  weak  projectile-target  interaction.  A  transient  regime  exists  at  intermediate 
velocities.  The  transition  regime  is  of  particular  interest  as  a  local  minimum  in  the  residual 
velocity  is  present.  In  order  to  capture  this  response  in  the  Lambert-Jonas  fonnulism,  a 
two-stage  expanded  Lambert-Jonas  model  is  defined 


V 


0,  V„p  <  VM 

■a/C-Cfi1.  VM  <¥„„<¥, 

’  V^>V« 


(0) 


with  Virans  the  transition  velocity  between  the  two  limiting  regimes,  Vy  is  a  virtual  ballistic 
velocity  fitted  on  high  impact  velocity  regime,  and  a\,  <22,  Pi  and  /Mhc  characterizing 
parameters.  Figure  22  depicts  the  residual  velocity  date  fitted  to  the  two-stage  model  with  fit 
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parameters  a\  =  0.308,  a2  =  0.425,  p\  =  1.772,  p2  =  1.528,  Vy  =  46  m/s,  and  Vtrans=  56  m/s, 
respectively.  The  coefficient  of  detennination  for  low  impact  velocity  is  0.9498  while  for 
high  impact  velocity  is  0.9883,  which  indicates  a  much  better  fitting  with  two-stage 
Lambert-Jonas  model  than  the  basic  model  for  the  high  impact  velocity  regime. 

To  further  explore  the  finding  of  the  regime  with  a  local  minimum  in  residual  velocity, 
computations  were  performed  with  a  lower  and  higher  value  of  p,  0.15  and  0.25  respectively. 
Figs. 15  and  16  depict  the  projectile  residual  velocities,  as  well  as  generalized  model  fits.  The 
fit  parameters  are  shown  Table  1.  In  fitting  the  parameters,  p\  and  pi  are  constrained  to  be 
greater  than  one,  as  demonstrated  in  [30].  Comparison  of  the  results  in  Figs.  14-16  shows  that 
the  ballistic  velocities  and  transition  velocities  increase  with  a  increasing  of  pL.  The 
charactering  parameters  a\,  pi  and  Vy  increase  with  p  while  ai  and  p\  decrease  monotonically. 
The  two  regimes  in  the  residual  velocities  are  present  independent  of  the  friction  condition, 
but  importance  of  the  transition  regime  with  minimal  velocity  diminishes  or  vanishes  as  the 
coefficient  of  friction  is  increased. 


Figure  22.  Predicted  residual  velocities  and  fit  by  two-stage  Lambert-Jonas  model  (ju=  0.20). 
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Figure  23.  Predicted  residual  velocities  and  fit  by  two-stage  Lambert-Jonas  model  (jl=  0.15). 


Figure  24.  Predicted  residual  velocities  and  fit  by  two-stage  Lambert-Jonas  model  (ji=  0.25). 
Table  1.  Parameters  of  the  two-stage  Lambert-Jonas  model 
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a\ 

Cl2 

Pi 

P2 

C, /[m/s] 

fV[m/s] 

V frans 

[m/s] 

R2 

p=  0.15 

0.240 

0.484 

2.230 

1.266 

28 

44 

55 

0.9083/0.9592 

ju=  0.20 

0.308 

0.425 

1.772 

1.528 

34 

46 

56 

0.9498/0.9883 

ju=  0.25 

0.325 

0.370 

1.721 

2.016 

39 

52 

65 

0.9914/0.9923 

4.3.  Fragmentation 

Fragmentation  commonly  occurs  for  brittle  materials  subject  to  impact  loading  [37-39].  To 
compare  the  impact  resistance  of  TIMs  made  by  brittle  materials  with  and  without 
fragmenting,  a  cohesive  zone  model  (CZM)  was  further  developed.  Motivated  by  the  TIM 
material  used  in  experiments  [15],  in  which  the  tetrahedra  were  made  with  additive 
manufacturing,  the  zero  thickness  cohesive  zone  (CZ)  layers  were  created  parallel  to  the 
mid-plane  of  TIM  assembly,  i.e.,  perpendicular  to  the  impact  direction,  which  was  meshed 
with  cohesive  element  COH3D8.  In  each  tetrahedron,  equally  spaced  eight  CZ  layers  were 
inserted.  Therefore,  each  tetrahedron  contains  10  solid  layers  and  8  CZ  layers  (Fig.  4).  A 
bilinear  traction-separation  law  was  used  to  describe  the  constitutive  response  of  CZM,  with 
maximum  strength  <Jmax=  El  100,  failure  separation  df=sol 1 000  and  damage  initiation 
separation  <^i=  <5f  / 1 0 .  Since  ABAQUS/Explicit  method  is  used  to  solve  the  dynamic 
problem,  the  density  of  the  CZ  layer  needs  to  be  defined.  In  the  simulations  for  TIMs  with 
CZ  layers,  the  density  of  the  CZ  layer  is  pcoh  =1.0  kg/m3  which  does  not  increase  the  inertia 
of  the  structure  significantly,  as  well  as  the  computational  difficulty  [34].  Aside  from  of  the 
CZ  properties,  all  the  other  properties  and  conditions  are  the  same  with  that  used  in  the 
impacting  simulations  of  TIMs  with  isotropic  tetrahedra. 

Four  impact  velocity  cases  ( Vimp  =  34,  45,  55  and  80  m/s)  were  perfonned  with  the  model 
enhanced  by  cohesive  zone  elements.  Table  2  lists  the  tenninal  velocities  for  isotropic  TIM 
and  CZ  enhanced  TIM  models.  Figure  25  compares  the  projectile  velocities  with  Vimp=55  m/s, 
and  similar  comparison  can  be  obtained  for  other  impact  velocity  cases  and  Figure  26  shows 
the  fragmentation  of  CZ  enhanced  TIMs  for  different  impact  velocities.  The  results  show 
that  the  CZ  enhanced  model  resists  impact  better  with  lower  projectile  residual  velocity,  but 
the  difference  is  not  significant  when  impact  velocity  is  high. 
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Table  2.  Terminal  velocities  of  TIMs  with  and  without  CZ  enhancement  [m/s] 


Cm/, =3  4  m/s 

Vimp=  45  m/s 

Vimp=55  m/s 

Vimp=  80  m/s 

Isotropic  TIM 

0 

8 

10 

23 

CZ  enhanced  TIM 

-2 

2 

8 

18 

Figure  25.  Projectile  velocity  for  isotropic  and  CZ  enhanced  TIM  models  ( Vimp=  55  m/s). 
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Figure  26.  Fragmentation  of  CZ  enhanced  TIMs  at  time  5  ms:  (a)  Vi,np  =  34  m/s  (b)  Vimp  = 

45  m/s  (c)  Vimp  =  55  m/s  (d)  Vimp  =  80  m/s. 

We  furthermore  explore  the  multiscale  characteristics  of  the  failure  response  of  TIMs  by 
expanding  the  considerations  on  the  cohesive  zone  model  approach  to  failure  in  the  context  of 
TIMs.  In  addition  to  the  parameter  set  considered  for  the  results  shown  in  the  initial  part  of 
this  section,  a  range  of  computations  was  conducted  in  which  the  cohesive  zone  model 
parameters  were  altered  systematically.  Table  3  summarizes  the  failure  parameters.  In  the 
model  series  the  internal  material  length  scale  and  the  material  strength  are  altered,  still  the 
material  stiffness  and  the  failure  energy  are  kept  constant.  Such  a  series  of  computations 
allows  for  the  study  of  the  influence  of  the  internal  material  length  scale  on  the  failure 
behavior.  Figure  27  depicts  the  predicted  residual  velocities  of  the  impactor  (initial  velocity 
Vimp  =  80  m/s)  after  interaction  with  the  target. 


25 


Table  3:  Cohesive  Zone  Model  parameters 


Test  # 

G 

max 

4 

L  Is 

c  0 

CZ  Energy 

1 

El  50 

5  /2000 

0 

0.025 

s  E/100000 

0 

2 

El  100 

5  /1000 

0 

0.1 

5  E/100000 

0 

3 

E/200 

5  /500 

0 

0.4 

5  E/100000 

0 

4 

E/316.2 

5  /3 16.2 

0 

1 

5  E/100000 

0 

5 

E/400 

5  /250 

0 

1.6 

5  E/100000 

0 

6 

E/447.2 

5  /223.614 

0 

2 

5  E/100000 

0 

7 

E/500 

5  /200 

0 

2.5 

5  E/100000 

0 

8 

E/1000 

5  /100 

0 

10 

5  E/100000 

0 

9 

E/2000 

5  /50 

0 

40 

5  E/100000 

0 

(Material  Length  scale  Lc)/(Unit  element  size  s„) 

Figure  27:  Predicted  residual  velocity  of  the  penetrator  in  dependence  of  the  material  failure 

characteristic  length  scale. 
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The  results  of  Fig.  27  indicate  a  size  effect  on  the  residual  velocity  in  dependence  of  the 
material  characteristic  length  scale.  For  cases  where  the  material  length  scale  is  large 
compared  to  the  size  of  the  unit  tetrahedron,  the  residual  velocity  was  found  to  be  nearly 
constant  and  independent  of  the  material  length  scale. 

5.  Discussion 

The  simulations  of  the  TIM  response  under  constant  loading  speed,  Figure  14  to  16,  provide 
insight  into  the  relevance  of  the  material  parameters  in  the  TIM  system.  The  responses  can  be 
divided  into  three  categories  according  to  loading  speeds:  quasi-static  response  at  low  speeds, 
combination  of  quasi-static  and  dynamic  response,  and  dominant  dynamic.  For  quasi-static 
case,  a  quasi-ductile  response  can  be  observed,  which  agrees  with  the  experimental 
observation  [15].  For  moderate  indentation  speeds,  the  force-deflection  response  can  be 
considered  as  the  combination  of  quasi-static  response  and  unconfined  TIM  under  same 
indention  speed.  As  shown  in  Appendix,  the  transition  speed  is  detennined  by  the  density  of 
TIM  material,  the  normal  contact  behavior,  while  coefficient  of  friction  does  not  affect  the 
inertia  responses  significantly. 

The  dynamic  characteristics  of  TIM  are  further  investigated  with  projectile  impact  loading 
conditions.  The  projectile  residual  velocity  is  used  as  the  parameter  to  illustrate  TIM’s 
resistibility  to  impact,  with  Lambert-Jonas  models  being  employed  to  describe  the  numerical 
results.  To  have  a  more  accurate  description  of  the  residual  velocity  results,  an  extended 
two-stage  Lambert-Jonas  model  is  developed.  For  conventional  plate,  Borvik  et  al.  [40]  fitted 
the  Lambert-Jonas  formula  with  a  =  0.76  and  p  =  2.36based  on  experimental  data  for  ballistic 
penetration  of  steel  plates.  Chi  et  al.  [41]  applied  the  model  to  a  bi-layer  metal  system  and 
obtained  ranges  of  the  parameters  for  different  thicknesses  of  metal  plates:  0.93  <a<  0.97  and 
1.72  <  p  <  2.14.  Other  experimental  and  numerical  results  [42-49]  also  indicate  that  for 
homogeneous  plates  and  composite  plates,  a  is  close  to  1.0  and  p  is  close  to  2.0.  Comparing 
to  the  parameters  in  the  Lambert-Jonas  model  for  TIM  with  conventional  plates,  the 
parameter  p  is  close  to  that  for  conventional  plates  while  a  is  significantly  smaller.  According 
to  Lambert-Jonas  formulism,  a  smaller  a  indicates  lower  residual  velocity,  therefore,  more 
kinetic  energy  absorption. 

Residual  velocity  comparison  was  performed  between  TIMs  with  isotropic  tetrahedra  or  CZ 
enhanced  tetrahedra.  When  the  impact  velocity  is  close  to  ballistic  velocity  of  TIM  with 
isotropic  tetrahedra,  the  residual  velocities  for  TIM  with  CZ  layers  are  smaller.  The 
increasing  of  TIM  stiffness  is  due  to  the  fragmentation.  At  the  end  of  impact  events,  energy 
distribution  of  the  entire  system  show  that,  the  kinetic  energy  of  the  tetrahedra  and  fragments 
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of  TIMs  with  CZ  layers  is  higher  than  TIMs  with  isotropic  tetrahedra.  These  results  can  be 
confirmed  with  the  result  shown  in  Figure  26.  For  example,  when  the  impact  velocity  is  34 
m/s,  every  tetrahedron  is  almost  motionless  for  TIMs  without  CZ  layers  (Figure  20(a))  at  the 
end  of  impact;  while  for  TIMs  with  CZ  layers  (Figure  26(a)),  fragments  fly  with  certain 
kinetic  energy.  Therefore,  because  of  fragmentation,  the  implementation  of  CZ  layers  into 
TIM  actually  increases  TIM's  dynamic  stiffness.  The  increment  in  dynamic  stiffness  is 
remarkable  when  impact  velocities  are  close  to  ballistic  velocity.  For  high  impact  velocities 
( V i,np  =  55  and  80  m/s),  the  difference  of  kinetic  energies  of  tetrahedra/fragments  for  both 
models  is  not  significant,  resulting  in  similar  residual  velocities  as  shown  in  Table.  2. 
Figure  26  also  illustrates  that  damage  can  be  localized  for  TIMs:  only  the  tetrahedra  directly 
impacted  and  their  neighbors  are  damaged,  while  the  tetrahedra  away  from  the  impact  center 
keep  their  integrities.  A  critical  internal  length  scale  appears  to  exist.  For  smaller  internal 
material  length  scales  that  those  of  the  critical  dimension,  the  residual  velocity  is  found  to  be 
dependent  on  the  material  length  scale  and  is  predicted  to  increase  with  a  decrease  in  material 
length  scale.  This  finding  appears  to  parallel  the  size  effects  found  in  quasi-brittle  structures 
as  widely  discussed  in  the  structural  mechanics  literature. 


6.  Conclusions 

The  present  work  numerically  investigated  the  dynamic  responses  of  TIM  under  constant 
velocity  and  projectile  impact  loading  with  finite  element  method.  The  FEA  model  is 
calibrated  with  published  experimental  data.  The  resistances  of  TIM  under  constant  velocity 
and  projectile  impact  were  examined.  Numerical  results  indicate  that:  (i)  TIM  responses  to 
constant  velocity  loading  can  be  categorized  into  three  groups  according  to  the  indention 
speeds:  quasi-static  response,  combination  of  dynamic  and  quasi-static  response,  and  pure 
dynamic  response.  At  low  velocities,  the  response  of  TIMs  emerges  as  the  sum  of  the  inertia 
(unconfined  TIM)  and  the  quasi-static  response.  For  higher  velocity  cases,  the  indenter 
imparts  multiple  inertia  responses.  (ii)The  terminal  velocities  of  projectile  for  impact 
velocities  form  1  m/s  to  80  m/s  are  obtained,  Investigation  of  the  residual  velocities  of 
projectile  after  penetrating  shows  TIMs  follow  two-stage  generalized  analytical  models.  The 
parameters  in  extended  Lambert-Jonas  models  for  TIM  are  fitted  on  obtained  numerical 
results,  with  parameter  p  is  close  to  conventional  solid  metallic  plates  for  high  velocity 
regime  but  is  smaller  for  low  velocity  regime,  while  parameter  a  is  much  smaller  for  both 
stages.  For  given  p  in  Lambert-Jonas  formulism,  a  smaller  a  indicates  lower  residual 
velocity,  hence  higher  kinetic  energy  absorption,  (iii)  Comparison  between  TIMs  with  and 
without  CZ  layers  show  that  TIMs  with  CZ  layers  can  resist  projectile  better  when  the  impact 
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velocity  is  close  to  ballistic  velocity.  When  the  impact  velocity  is  higher,  the  difference  of  the 
two  TIM  models  in  resisting  impact  is  not  significant.  Here  we  demonstrate  a  size  effect  on 
failure  for  the  topologically  interlocked  materials  and  provide  a  path  such  that  these  material 
systems  could  be  exploited  to  create  impact  resistant  systems. 


7.  Collaboration 

During  the  project  period  the  PI  has  conducted  extensive  coordination  and  communication 
with  research  staff  at  ARL,  E.  Habtour  and  J.  Riddick.  The  PI  has  provided  ARL  staff  the 
Python  scripting  code,  stl  files  for  3D  printing  of  TIM  assemblies  in  house  at  ARL,  and  has 
delivered  input  files  for  ABAQUS  computations  to  ARL.  The  PI  visited  ARL  on  Sep.  26, 
2013  and  presented  a  seminar  talk  summarizing  the  results  of  this  work. 


8.  Publications 

The  following  publications  have  resulted  from  this  effort: 

Journal  Publications 

Y.-Z.  Peng,  E.  Habtour,  J.  Riddick,  T.  Siegmund,  “Impact  mechanics  of  topologically 
interlocked  materials,”  International  Journal  of  Impact  Engineering,  to  be  submitted. 

ARL-Reports 

Y.-Z.  Feng,  T.  Siegmund,  “Python  scripting  in  ABAQUS  to  create  finite  element  models  and 
3D  print  physical  models  for  granular  crystals,”  submitted  to  ARL  as  Technical  Report. 

Conferences  Presentations: 

Y.-Z.  Feng,  E.  Habtour,  J.  Riddick,  T.  Siegmund,  “Damage  analysis  of  a  hybrid  energy 
absorption  layer  based  on  the  principle  of  topologically  interlocking  materials,”  Proceedings 
of  the  2013  International  Mechanical  Engineering  Congress  and  Exposition,  The  American 
Society  of  Mechanical  Engineers.  San  Diego,  CA,  Nov.  2013;  paper  IMECE20 13-66426. 
Y.-Z.  Feng,  T.  Siegmund,  “Dynamic  and  failure  analysis  of  topologically  interlocked 
materials  subjected  to  impact,”  EMI-2013,  2013  Engineering  Mechanics  Institute  Conference 
of  ASCE,  August  4-7,  2013,  Northwestern  University,  Evanston,  IL,  USA,  Abstract  #371. 
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Appendix.  Determination  of  the  TIM  Inertia 

To  investigate  the  determination  of  the  ’inertia'  of  TIMs,  unconfined  TIM  models  were 
employed  for  different  tetrahedron  material  densities,  COFs  and  contact  stiffness.  Loading 
speed  is  Vo  =  4  m/s  for  all  the  unconfined  cases.  Figure  A-l  shows  the  predicted  peak  forces 
for  different  material  densities.  The  material  densities  are  normalized  with  the  density  of 
engineering  ABS  (po  =  950  kg/m3).  Approximate  linear  dependency  of  peaks  with  the  density 
can  be  observed  from  Figure  A-l.  When  p  =  0.05p0,  the  peak  of  the  reaction  forces  is  about 
20  N,  while  when  p  =  2.0 po,  the  peak  of  the  reaction  forces  is  about  177  N.  Figure  A-2  shows 
the  predicted  peak  forces  for  normalized  COFs  range  from  0.05  to  2.0.  The  COFs  are 
normalized  with  calibrated  COF  (po  =  0.2).  With  40  times  increase  (2/0.05)  in  COF,  the  peak 

value  increases  only  about  33%.  Figure  A-3  depicts  the  predicted  peak  forces  for  different 

* 

contact  stiffness.  The  contact  stiffness  is  normalized  with  calibrated  value  ( K  o  =  0.38).  The 
peak  force  is  about  39  N  when  K*=  0.05A*o,  and  reaches  154  N  when  K*=  2.0K*o,  a  total  4 
time  increase  over  the  range  of  contact  stiffness  from  0.05A'o  to  2.0  Kq. 


Figure  A-l.  Predicted  peak  reaction  forces  of  unconfined  TIM  for  various  values  material 

densities. 
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Figure  A-2.  Predicted  peak  reaction  forces  of  unconfined  TIM  for  various  values  of 

coefficients  of  friction. 


Figure  A-3.  Predicted  peak  reaction  forces  of  unconfined  TIM  for  various  contact  stiffness. 
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List  of  Symbols  and  Abbreviations 

a:  Parameter  in  Lambert- Jonas  model 

a\,  ay  Parameters  in  extended  Lambert- Jonas  models 

Cd.  Dilatational  wave  speed 

D :  Projectile  diameter 

E:  Young's  modulus  of  tetrahedron  material 

[I drop-  Drop  height  in  drop  tower  testing 

K :  Contact  stiffness 

K*:  Normalizedcontact  stiffness 

L\  Mid-plane  dimension  of  topologically  interlocked  material  assembly 
mo'.  Mass  of  a  single  tetrahedron 
M\  Projectile  mass 

Mdrop :  Drop  mass  in  drop  tower  testing 

N:  Number  of  tetrahedra  in  one  row/column  of  topologically  interlocked  material  assembly 
p\  Parameters  in  Lambert-Jonas  models 

p\  and  py  Parameters  in  extended  two-stage  Lambert-Jonas  models 

.so:  Edge  length  of  regular  tetrahedron 

t:  Thickness  of  topologically  interlocked  material  assembly 

Vo:  Constant  velocity 

Vbf.  Ballistic  velocity 

Vimp:  Impact  velocity 

Vres:  Residual  velocity 

Vv :  Fitted  virtual  ballistic  velocity  of  TIM 

<5i),  8f.  Damage  initiation  separation  and  failure  separation  in  CZM 
v:  Poison's  ratio  of  tetrahedron  material 
p:  Coefficient  of  friction 

p0,p:  Densities  of  tetrahedron  and  projectile  respectively 
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pcoh.  Density  of  cohesive  layers 

Omax’  Maximum  strength  of  CZM 

COF:  Coefficient  of  friction 

CZ,  CZM:  Cohesive  zone  and  cohesive  zone  model 

TIM:  Topologically  interlocked  material 
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